% Changing C and gL from pyr to FS values made almost no difference in the
% gain of the cell.
function [spikes Gampa first_spike_rate]=adaptexp(stim_conductance_exc, stim_conductance_inh,neuron,g_true_I_false)

global v I

tauVT=1; %default;
DeltaVT=0; % default;
if isstruct(neuron)
    Params = neuron;
        % unbundle structure
    C = Params.C;
    gL = Params.gL;
    EL = Params.EL;
    VT = Params.VT;
    DeltaT = Params.Delta;
    Vcut=VT+Params.n*DeltaT
    a = Params.a;
    tauw = Params.TauW;
    b = Params.b;
    Vr = Params.Er;
    if isfield(Params,'tauVT')
        tauVT = Params.tauVT;
        DeltaVT=Params.DeltaVT;
    end
    
else

    switch neuron
        % pyr
        case 1
        C=104; % *pF
        gL=4.3; %*nS
        taum=C/gL; 
        EL=-65; %*mV # Same as changing I
        VT=-52;%-52; %*mV  -56 gives a chattering cell
        DeltaT=.8; %*mV
        Vcut=VT+5*DeltaT;
        b=65; %*nA
        Vr=-53; %*mV
        tauw=88;
        a=-.8; %-.8;

        case 2

        % FS cell 
        C=59; % *pF
        gL=2.9; %*nS
        taum=C/gL;
        EL=-62; %*mV # Same as changing I
        VT=-42; %*mV
        %VT=-32 %*mV
        DeltaT=3; %*mV
        Vcut=VT+5*DeltaT;
        b=61; %*nA
        Vr=-54; %*mV
        tauw=16;
        a=1.8;


    end
end

VTa_decay = exp(-1/tauVT);

if ~exist('g_true_I_false')
    g_true_I_false = true;
end

n=max(length(stim_conductance_inh), length(stim_conductance_exc))
   
if n==1

    T=1000;

    dt_ms=1;
    n=round(T/dt_ms);
else
    dt_ms=1;
end

spikes = 0;

v=zeros(1,n);
u=zeros(1,n);
v(1)=EL;
u(1)=0;


stp_tau = exp(-dt_ms/150);
stp = 0*stp_tau.^(0:n-1)+1;

si= stim_conductance_inh.*[ones(1,n)];
se= stim_conductance_exc.*[ones(1,n)];

% model conductances.
tau_ampa  = exp(-dt_ms/5);
tau_nmda  = exp(-dt_ms/150);
tau_gabaa = exp(-dt_ms/6);
tau_gabab = exp(-dt_ms/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70; % Niraj: -65;
E_GABAb = -90;
Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;

Ggabaa=zeros(1,n);
Ggabab=zeros(1,n);
Gampa=zeros(1,n);
Gnmda=zeros(1,n);
I=zeros(1,n);

gain_ampa=1;
gain_nmda=.1;
gain_gabaa=1;
gain_gabab=.1;
% end model conductances.

first_spike_rate = 0;
VT_adaption = 0;
for i = 2:n
    
    % conductances from synapses.
    Ggabaa(i) = si(i);
    %Ggabab(i) = Ggabab(i-1)*tau_gabab + (1-tau_gabab)*gain_gabab*si(i);
    
    Gampa(i) = se(i);
    %Gnmda(i) = Gnmda(i-1)*tau_nmda + (1-tau_nmda)*gain_nmda*se(i);

    xx = (-80-v(i))/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
     
    if g_true_I_false
        g= ( Gampa(i)       +NMDAgate.*Gnmda(i)        + Ggabaa(i)         + Ggabab(i));
        E= ( Gampa(i)*E_AMPA+NMDAgate.*Gnmda(i)*E_NMDA + Ggabaa(i)*E_GABAa + Ggabab(i)*E_GABAb);
        I(i) = v(i-1).*g-E;
        %I(i) = (v(i-1)-E_AMPA)*Gampa(i) + (v(i-1)-E_NMDA)*NMDAgate*Gnmda(i)+ ...
        % (v(i-1)-E_GABAa)*Ggabaa(i)+(v(i-1)-E_GABAb)*Ggabab(i);
    else
        I(i) =se(i);
    end
    
    dvm=(gL*(EL-v(i-1))+gL*DeltaT*exp((v(i-1)-(VT+VT_adaption))/DeltaT)-I(i)-u(i-1))/C; % volt
    v(i) = v(i-1) + dt_ms*dvm;
    if v(i) < -90
        v(i) = -90;
    end
    if v(i) > 51
        v(i) = 51;
    end
    dw=(a*(v(i)-EL)-u(i-1))/tauw; % amp
    u(i)=u(i-1) + dt_ms*dw;
    
    dVTa = -VT_adaption/tauVT;
    VT_adaption = VT_adaption + dt_ms*dVTa;
    
    if( v(i) > Vcut )
        v(i-1)=0;
        v(i) = Vr;
        u(i) = u(i)+b;
        spikes = spikes + 1;
        if spikes == 1
            first_spike_ms = i;
        elseif spikes == 2
            isi1=i-first_spike_ms;
            first_spike_rate = 1000/isi1;
        end
        
        VT_adaption = VT_adaption + DeltaVT; 
    end

   
end

figure(3)
subplot(3,1,1)
plot(dt_ms*(1:n),[v' u']);
subplot(3,1,2)
plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
title('conductances');
subplot(3,1,3);
plot(dt_ms*(1:n),I);
title('Isyn (pA)');
drawnow
